#!/usr/bin/python3.9
# -*- coding: utf-8 -*-
# @Time    : 2021/10/10 17:16
# @Author  : YHSimon
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.io import loadmat
from scipy.optimize import minimize


def load_data(path):
    data = loadmat(path)
    X = data['X']
    y = data['y']
    return X, y


def plot_an_image(X):
    # 随机打印一个数字
    pick_one = np.random.randint(0, 5000)
    image = X[pick_one, :]
    fig, ax = plt.subplots(figsize=(1, 1))
    ax.matshow(image.reshape((20, 20)), cmap='gray_r')
    plt.xticks([])  # 去除梯度，美观
    plt.yticks([])
    plt.show()
    print('this should be {}'.format(y[pick_one]))


def plot_100_image(X):
    # 随机画100个数字
    sample_idx = np.random.choice(np.arange(X.shape[0]), 100)  # 随机选100个样本
    sample_images = X[sample_idx, :]  # (100,400)   20×20的像素网格被展开成一个400维的向量
    # sharex和sharey表示坐标轴的属性是否相同  nrows和ncols表示将画布分割成几行几列
    fig, ax_array = plt.subplots(nrows=10, ncols=10, sharey=True, sharex=True, figsize=(8, 8))
    for row in range(10):
        for col in range(10):
            ax_array[row, col].matshow(sample_images[10 * row + col].reshape((20, 20)), cmap='gray_r')
    plt.xticks([])
    plt.yticks([])
    plt.show()


def sigmoid(z):
    return 1 / (1 + np.exp(-z))


def regularized_cost(theta, X, y, l):
    """ 向量化代价函数
    不惩罚θ_0
    :param theta:
    :param X: 特征矩阵 （m,n+1） # 插入了x0=1
    :param y: 目标向量  (m,)
    :param l: 正则化λ常熟
    """
    thetaReg = theta[1:]
    first = (-y * np.log(sigmoid(X @ theta))) + (y - 1) * np.log(1 - sigmoid(X @ theta))
    reg = (thetaReg @ thetaReg) * l / (2 * len(X))
    return np.mean(first) + reg


def regularized_gradient(theta, X, y, l):
    """
    不惩罚θ_0
    :param theta:
    :param X:
    :param y:
    :param l: 正则化常数λ
    """
    thetaReg = theta[1:]
    first = (1 / len(X)) * X.T @ (sigmoid(X @ theta) - y)
    # 这里人为插入一组0，使得theta_0不惩罚，方便计算
    reg = np.concatenate([np.array([0]), (l / len(X)) * thetaReg])
    return first + reg


"""
    这部分我们将实现一对多分类通过训练多个正则化logistic回归分类器，每个对应数据集中K类中的一个。对于这个任务，我们有10个可能的类，
    并且由于logistic回归只能一次在2个类之间进行分类，每个分类器在“类别 i”和“不是 i”之间决定。 我们将把分类器训练包含在一个函数中，
    该函数计算10个分类器中的每个分类器的最终权重，并将权重返回shape为(k, (n+1))数组，其中 n 是参数数量。
    
    这里需要注意的几点：首先，我们为X添加了一列常数项 1 ，以计算截距项（常数项）。 其次，我们将y从类标签转换为每个分类器的二进制值（要么是类i，要么不是类i）。 
    最后，我们使用SciPy的较新优化API来最小化每个分类器的代价函数。 如果指定的话，API将采用目标函数，初始参数集，优化方法和jacobian（渐变）函数。 
    然后将优化程序找到的参数分配给参数数组。

"""


def one_vs_all(X, y, l, K):
    """
    :param X: feature matrix, (m, n+1) # with incercept x0=1
    :param y: target vector, (m, )
    :param l: lambda constant for regularization
    :param k: numbel of labels
    """
    all_theta = np.zeros((K, X.shape[1]))  # (10,401)
    for i in range(1, K + 1):
        theta = np.zeros(X.shape[1])
        y_i = np.array([1 if label == i else 0 for label in y])
        ret = minimize(fun=regularized_cost, x0=theta, args=(X, y_i, l), method='TNC',
                       jac=regularized_gradient, options={'disp': True})
        all_theta[i - 1, :] = ret.x
    return all_theta


def predict_all(X, all_theta):
    # compute the class probability for each class on each training instance
    h = sigmoid(X @ all_theta.T)  # 注意的这里的all_theta需要转置
    # create array of the index with the maximum probability
    # Returns the indices of the maximum values along an axis.
    h_argmax = np.argmax(h, axis=1)
    # because our array was zero-indexed we need to add one for the true label prediction
    h_argmax = h_argmax + 1

    return h_argmax  # 返回的h_argmax是一个array，包含5000个样本对应的预测值。


X, y = load_data('ex3data1.mat')
print(np.unique(y))  # 看下有几类标签
# [ 1  2  3  4  5  6  7  8  9 10]
print(X.shape, y.shape)  # ((5000, 400), (5000, 1))
plot_an_image(X)
plot_100_image(X)

if __name__ == '__main__':
    raw_X, raw_y = load_data('ex3data1.mat')
    X = np.insert(raw_X, 0, 1, axis=1)  # (5000, 401)
    y = raw_y.flatten()  # 这里消除了一个维度，方便后面的计算 or .reshape(-1) （5000，）
    # y = raw_y.reshape(-1)

    all_theta = one_vs_all(X, y, 1, 10)  # 每一行是一个分类器的一组参数
    print(all_theta)
    y_pred = predict_all(X, all_theta)
    accuracy = np.mean(y_pred == y)
    print('accuracy = {0}%'.format(accuracy * 100))  # accuracy = 94.46%
